set more off
**read in dataset
use data/bat_dataset, clear
/* variable definitions

fips: county identifier
year: year of the observation
dryrent: county rental rate for non-irrigated land
agland: number of acres in cropland
WNS: indicator equal to 1 if a county has had an outbeak of WNS at or before year
neighbor_outbreak_ind: indicator equal to 1 if at least 1 neighbor had an outbreak of WNS in or before year
neighbor_outbreak_count: count of the number of neighbors that had an outbreak of WNS in or before year
q?precip: precip in quarter ? of year
q?mintemp: average daily min temperature in quarter ? of year
pop: human population
st: state
rg: region
lncyield: natural log of corn yields
lnwyield: natural log of wheat yields
lnsyield: natural log of soy yields
exp_land: input expenditure per acre
highqual: proportion of cropland that is high quality
medqual: proportion of cropland that is medium quality
lowqual: proportion of cropland that is low quality
qual1-qual4: proportion of cropland in LCC class 1-4
quallow: proportion of cropland in LCC classes>4
lnagland: natural log of the number of acres in cropland
lndryrent: natural log of dryrent
corn_yield: corn yield (bu/acre)
wheat_yield: wheat yield (bu/acre)
soy_yield: soy yield (bu/acre)
wn: indicator equal to 1 if a county ever had a WNS outbreak
event: time period of the analysis when WNS first equals 0 (used for Figure 3)
event2: number of years since first outbreak (0 in year of first outbreak)  (used for Figure 3)
*/
xtset fips year	
	********establish samples from main models:
set more off
**rent model sample
xtreg dryrent WNS neighbor_outbreak_ind q?precip q?mintemp pop i.year#i.st if year >=1980 & rg !=6 & rg !=4 , fe cluster(fips)
	gen sample = e(sample)
	bys fips: egen insample = max(sample)
	gen WNS_sample = WNS if e(sample)
	gen WNS_n_count_sample = neighbor_outbreak_count if e(sample)
	gen WNS_n_n_sample = neighbor_outbreak_ind  if e(sample)
**agland model sample
xtreg lnagland WNS neighbor_outbreak_count  q?precip q?mintemp pop  i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample==1  , fe cluster(fips)
	gen agland_sample = agland if e(sample)	
	
***Table 2***********	
xtreg dryrent WNS neighbor_outbreak_count  q?precip q?mintemp pop highqual medqual   i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1, fe cluster(fips)
	gen dryrent_sample = dryrent if e(sample)
	outreg2 using output\Table2,  replace keep (WNS neighbor_outbreak_count ) tex
xtreg dryrent WNS neighbor_outbreak_ind  q?precip q?mintemp pop highqual medqual   i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1, fe cluster(fips)
	outreg2 using output\Table2,   keep (WNS neighbor_outbreak_ind ) tex
xtreg dryrent WNS neighbor_outbreak_count  i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1, fe cluster(fips)
	outreg2 using output\Table2,   keep (WNS neighbor_outbreak_count ) tex
xtreg lndryrent WNS neighbor_outbreak_count  q?precip q?mintemp pop highqual medqual   i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1, fe cluster(fips)
	outreg2 using output\Table2,   keep (WNS neighbor_outbreak_count ) tex
xtreg dryrent WNS   q?precip q?mintemp pop highqual medqual   i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1, fe cluster(fips)
	outreg2 using output\Table2,  keep (WNS neighbor_outbreak_count ) tex
*********************************************************************************	
**********Table 3 (not all directly output using outreg2)************
preserve
	xtreg dryrent WNS neighbor_outbreak_count q?precip q?mintemp pop i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1 ,  fe  cluster(fips)
			scalar beta = _b[WNS]
			scalar delta = _b[neighbor_outbreak_count]
		gen newvar = WNS	+delta/beta * neighbor_outbreak_count
	xtreg agland  newvar q?precip q?mintemp pop i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1,  fe cluster(fips) 
		scalar gamma_main = _b[newvar]
		scalar alpha_main = _b[_cons]
		lincom _b[newvar]*delta/beta
		scalar theta_main = r(estimate)	

	xtreg agland  newvar  i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1,  fe cluster(fips) 
		scalar gamma_noc = _b[newvar]
		scalar alpha_noc = _b[_cons]
		lincom _b[newvar]*delta/beta
		scalar theta_noc = r(estimate)	
	matrix test = (beta, delta, gamma_main, theta_main, alpha_main, gamma_noc, theta_noc, alpha_noc )
matrix list test
	**with bootstrap***
	capture program drop myboot
		program define myboot, rclass
		preserve		
			bsample, cluster(fips)
			set more off
			xtreg dryrent WNS neighbor_outbreak_count q?precip q?mintemp pop i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1 ,  fe  cluster(fips)
				return scalar beta = _b[WNS]
				return scalar delta = _b[neighbor_outbreak_count]
			capture drop newvar
			gen newvar = WNS	+delta/beta * neighbor_outbreak_count

		xtreg agland  newvar q?precip q?mintemp pop i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1,  fe cluster(fips) 
			return scalar gamma_main = _b[newvar]
			return scalar alpha_main = _b[_cons]
			lincom _b[newvar]*delta/beta
			return scalar theta_main = r(estimate)
			
		xtreg agland  newvar  i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1,  fe cluster(fips) 
			return scalar gamma_noc = _b[newvar]
			return scalar alpha_noc = _b[_cons]
			lincom _b[newvar]*delta/beta
			return scalar theta_noc = r(estimate)	
		restore
	end
	
	simulate beta = r(beta) delta= r(delta)   gamma_main = r(gamma_main) theta_main = r(theta_main) alpha_main = r(alpha_main) ///
	gamma_noc = r(gamma_noc) theta_noc = r(theta_noc) alpha_noc = r(alpha_main) , reps(200) seed(123456789): myboot
 bstat, stat(test) n(257) title("Bootstrap Results--Samples Drawn by County")

restore

quietly xtreg agland WNS  neighbor_outbreak_count   q?precip q?mintemp pop   i.year#i.st if year >=1980 & rg !=6 & rg !=4  & insample==1  , fe cluster(fips)
	outreg2 using output\Table3,   keep (WNS neighbor_outbreak_count  ) tex replace
quietly xtreg lnagland WNS  neighbor_outbreak_count   q?precip q?mintemp pop   i.year#i.st if year >=1980 & rg !=6 & rg !=4  & insample==1  , fe cluster(fips)
	outreg2 using output\Table3,   keep (WNS neighbor_outbreak_count  ) tex
quietly xtreg agland WNS    q?precip q?mintemp pop   i.year#i.st if year >=1980 & rg !=6 & rg !=4  & insample==1  , fe cluster(fips)
	outreg2 using output\Table3,   keep (WNS neighbor_outbreak_count  ) tex
	
*******************************************************************
*********table 4**********************************
*yields
	**with neighbor outbreak
xtreg lncyield WNS neighbor_outbreak_count q?precip q?mintemp pop highqual medqual  i.year#i.rg if year >=1980 & rg !=6 & rg !=4  & insample==1 , fe cluster(fips)
	gen corny_sample = corn_yield if e(sample)
	outreg2 using output\table4, replace keep (WNS  neighbor_outbreak_count ) tex
xtreg lnwyield WNS  neighbor_outbreak_count q?precip q?mintemp pop highqual medqual i.year#i.rg if year >=1980 & rg !=6 & rg !=4  & insample==1  , fe cluster(fips)
	gen wheaty_sample = wheat_yield if e(sample)
	outreg2 using output\table4,  keep (WNS neighbor_outbreak_count ) tex
xtreg lnsyield WNS neighbor_outbreak_count  q?precip q?mintemp pop highqual medqual i.year#i.rg if year >=1980 & rg !=6 & rg !=4  & insample==1   , fe cluster(fips)
	gen soyy_sample = soy_yield if e(sample)
	outreg2 using output\table4,  keep (WNS  neighbor_outbreak_count ) tex
**excluding neighbor outbreak
xtreg lncyield WNS q?precip q?mintemp pop highqual medqual  i.year#i.rg if year >=1980 & rg !=6 & rg !=4  & insample==1 , fe cluster(fips)
	outreg2 using output\table4,  keep (WNS  neighbor_outbreak_count ) tex
	set more off
xtreg lnwyield WNS  q?precip q?mintemp pop highqual medqual i.year#i.rg if year >=1980 & rg !=6 & rg !=4  & insample==1  , fe cluster(fips)
	outreg2 using output\table4,  keep (WNS neighbor_outbreak_count ) tex
	set more off
xtreg lnsyield WNS  q?precip q?mintemp pop highqual medqual i.year#i.rg if year >=1980 & rg !=6 & rg !=4  & insample==1   , fe cluster(fips)
	outreg2 using output\table4,  keep (WNS  neighbor_outbreak_count ) tex

***************************************************************************888
***************table 5
******************************************************************************
********input costs
gen lnexp_land = ln(exp_land)
*including neighbor
xtreg lnexp_land WNS neighbor_outbreak_count q?precip q?mintemp pop  highqual medqual  i.year#i.st if year >=1990 & rg !=6 & rg !=4  & insample==1 , fe cluster(fips)
	gen exp_land_sample  = exp_land if e(sample)
	outreg2 using output\table5 ,replace keep(WNS neighbor_outbreak_count) tex
xtreg lnexp_land WNS neighbor_outbreak_count q?precip q?mintemp pop highqual medqual    i.year#i.rg if year >=1990 & rg !=6 & rg !=4  & insample==1   , fe cluster(fips)
	outreg2 using output\table5 ,keep(WNS neighbor_outbreak_count) tex
xtreg lnexp_land WNS neighbor_outbreak_count i.year#i.st  if year >=1990 & rg !=6 & rg !=4 & insample==1    , fe cluster(fips)
	outreg2 using output\table5 ,keep(WNS neighbor_outbreak_count) tex
*excluding neighbor
xtreg lnexp_land WNS q?precip q?mintemp pop  highqual medqual  i.year#i.st if year >=1990 & rg !=6 & rg !=4  & insample==1 , fe cluster(fips)
	outreg2 using output\table5 , keep(WNS neighbor_outbreak_count) tex
xtreg lnexp_land WNS q?precip q?mintemp pop highqual medqual    i.year#i.rg if year >=1990 & rg !=6 & rg !=4  & insample==1   , fe cluster(fips)
	outreg2 using output\table5 ,keep(WNS neighbor_outbreak_count) tex
xtreg lnexp_land WNS i.year#i.st  if year >=1990 & rg !=6 & rg !=4 & insample==1    , fe cluster(fips)
	outreg2 using output\table5 ,keep(WNS neighbor_outbreak_count) tex

*****************
*****************************
*****************************
*****appendix table c5 with limited data*************
**only 2008 and later
xtreg lncyield WNS  neighbor_outbreak_count q?precip q?mintemp pop highqual medqual  i.year#i.rg if year >=2008 & rg !=6 & rg !=4  & insample==1 , fe cluster(fips)
	outreg2 using output\appendix_mech,  keep (WNS neighbor_outbreak_count  ) tex replace
	set more off
xtreg lnwyield WNS neighbor_outbreak_count  q?precip q?mintemp pop highqual medqual i.year#i.rg if year >=2008 & rg !=6 & rg !=4  & insample==1  , fe cluster(fips)
	outreg2 using output\appendix_mech,  keep (WNS neighbor_outbreak_count ) tex
	set more off
xtreg lnsyield WNS  neighbor_outbreak_count q?precip q?mintemp pop highqual medqual i.year#i.rg if year >=2008 & rg !=6 & rg !=4  & insample==1   , fe cluster(fips)
	outreg2 using output\appendix_mech,  keep (WNS neighbor_outbreak_count  ) tex		
*only since 2008	
xtreg lnexp_land WNS neighbor_outbreak_count q?precip q?mintemp pop  highqual medqual  i.year#i.st if year >=2008 & rg !=6 & rg !=4  & insample==1 , fe cluster(fips)
	outreg2 using output\appendix_mech , keep(WNS neighbor_outbreak_count) tex
xtreg lnexp_land WNS neighbor_outbreak_count q?precip q?mintemp pop highqual medqual    i.year#i.rg if year >=2008 & rg !=6 & rg !=4  & insample==1   , fe cluster(fips)
	outreg2 using output\appendix_mech ,keep(WNS neighbor_outbreak_count) tex
xtreg lnexp_land WNS neighbor_outbreak_count i.year#i.st  if year >=2008 & rg !=6 & rg !=4 & insample==1    , fe cluster(fips)
	outreg2 using output\appendix_mech ,keep(WNS neighbor_outbreak_count) tex
**********************************************88**********8
************table 6
gen reg = 0
	replace reg = 1 if rg == 1 | rg ==2
	replace reg = 2 if rg == 3 | rg ==5 | rg ==7
	replace reg = 3 if rg == 8 | rg ==9
	set more off
	***rent by region
	scalar drop _all
xtreg dryrent WNS neighbor_outbreak_count q?precip q?mintemp pop highqual medqual i.year#i.st if year >=1980 & reg ==1, fe cluster(fips)
	outreg2 using output\table6,  keep(WNS  neighbor_outbreak_count ) replace tex
		scalar beta1 = _b[WNS]
		scalar delta1 = _b[neighbor_outbreak_count]
		capture drop newvar1
		gen newvar1 = WNS	+delta1/beta1 * neighbor_outbreak_count
	summ WNS if e(sample) & WNS ==1
xtreg dryrent WNS neighbor_outbreak_count   q?precip q?mintemp pop highqual medqual i.year#i.st if year >=1980 & reg ==2, fe cluster(fips)
	outreg2 using output\table6,  keep(WNS neighbor_outbreak_count  )  tex
		scalar beta2 = _b[WNS]
		scalar delta2 = _b[neighbor_outbreak_count]
		capture drop newvar2
		gen newvar2 = WNS	+delta2/beta2 * neighbor_outbreak_count
	summ WNS if e(sample) & WNS ==1
xtreg dryrent WNS neighbor_outbreak_count   q?precip q?mintemp pop highqual medqual i.year#i.st if year >=1980 & reg ==3, fe cluster(fips)
	outreg2 using output\table6,  keep(WNS  neighbor_outbreak_count )  tex
	summ WNS if e(sample) & WNS ==1
		scalar beta3 = _b[WNS]
		scalar delta3 = _b[neighbor_outbreak_count]
		capture drop newvar3
		gen newvar3 = WNS	+delta3/beta3 * neighbor_outbreak_count
		
	
**********************************************************
*******table 7*************************
	
	***land by region
	set more off
xtreg agland newvar1   q?precip q?mintemp pop i.year#i.st if year >=1980 & reg ==1 & insample ==1, fe cluster(fips)
	outreg2 using output\table7,  keep(WNS  neighbor_outbreak_count newvar1) replace tex
	summ WNS if e(sample) & WNS ==1
xtreg agland newvar2   q?precip q?mintemp pop i.year#i.st if year >=1980 & reg ==2, fe cluster(fips)
	outreg2 using output\table7,  keep(WNS neighbor_outbreak_count newvar2  )  tex
	summ WNS if e(sample) & WNS ==1
xtreg agland newvar3   q?precip q?mintemp pop i.year#i.st if year >=1980 & reg ==3, fe cluster(fips)
	outreg2 using output\table7,  keep(WNS neighbor_outbreak_count newvar3  )  tex
	summ WNS if e(sample) & WNS ==1		
	
	
preserve
	capture drop newvar*
	xtreg dryrent WNS neighbor_outbreak_count q?precip q?mintemp pop highqual medqual i.year#i.st if year >=1980 & reg ==1, fe cluster(fips)
			scalar beta1 = _b[WNS]
			scalar delta1 = _b[neighbor_outbreak_count]
		gen newvar1 = WNS	+delta1/beta1 * neighbor_outbreak_count
	xtreg agland newvar1   q?precip q?mintemp pop i.year#i.st if year >=1980 & reg ==1, fe cluster(fips)
		scalar gamma1 = _b[newvar1]
		scalar alpha1 = _b[_cons]
		lincom _b[newvar1]*delta1/beta1
		scalar theta1 = r(estimate)	
		
	xtreg dryrent WNS neighbor_outbreak_count q?precip q?mintemp pop highqual medqual i.year#i.st if year >=1980 & reg ==2, fe cluster(fips)
			scalar beta2 = _b[WNS]
			scalar delta2 = _b[neighbor_outbreak_count]
		gen newvar2 = WNS	+delta2/beta2 * neighbor_outbreak_count
	xtreg agland newvar2   q?precip q?mintemp pop i.year#i.st if year >=1980 & reg ==2, fe cluster(fips)
		scalar gamma2 = _b[newvar2]
		scalar alpha2 = _b[_cons]
		lincom _b[newvar2]*delta2/beta2
		scalar theta2 = r(estimate)
		
	xtreg dryrent WNS neighbor_outbreak_count q?precip q?mintemp pop highqual medqual i.year#i.st if year >=1980 & reg ==3, fe cluster(fips)
			scalar beta3 = _b[WNS]
			scalar delta3 = _b[neighbor_outbreak_count]
		gen newvar3 = WNS	+delta3/beta3 * neighbor_outbreak_count
	xtreg agland newvar3   q?precip q?mintemp pop i.year#i.st if year >=1980 & reg ==3, fe cluster(fips)
		scalar gamma3 = _b[newvar3]
		scalar alpha3 = _b[_cons]
		lincom _b[newvar3]*delta3/beta3
		scalar theta3 = r(estimate)
*restore
	
	matrix test = (beta1, delta1, gamma1, theta1, alpha1, beta2, delta2, gamma2, theta2, alpha2, beta3, delta3, gamma3, theta3, alpha3)
	matrix list test
*****change eps in 2 places	
	**with bootstrap***
	capture program drop myboot
		program define myboot, rclass
		preserve
		**need to adjust to sample clusters
		bsample, cluster(fips)
		set more off
		xtreg dryrent WNS neighbor_outbreak_count q?precip q?mintemp pop highqual medqual i.year#i.st if year >=1980 & reg ==1, fe cluster(fips)
			return scalar beta1 = _b[WNS]
			return scalar delta1 = _b[neighbor_outbreak_count]
		capture drop newvar1
		gen newvar1 = WNS	+delta1/beta1 * neighbor_outbreak_count

		xtreg agland newvar1   q?precip q?mintemp pop i.year#i.st if year >=1980 & reg ==1, fe cluster(fips)
			return scalar gamma1 = _b[newvar1]
			return scalar alpha1 = _b[_cons]
			lincom _b[newvar1]*delta1/beta1
			return scalar theta1 = r(estimate)
			
		xtreg dryrent WNS neighbor_outbreak_count q?precip q?mintemp pop highqual medqual i.year#i.st if year >=1980 & reg ==2, fe cluster(fips)
			return scalar beta2 = _b[WNS]
			return scalar delta2 = _b[neighbor_outbreak_count]
		capture drop newvar2
		gen newvar2 = WNS	+delta2/beta2 * neighbor_outbreak_count

		xtreg agland newvar2   q?precip q?mintemp pop i.year#i.st if year >=1980 & reg ==2, fe cluster(fips)
			return scalar gamma2 = _b[newvar2]
			return scalar alpha2 = _b[_cons]
			lincom _b[newvar2]*delta2/beta2
			return scalar theta2 = r(estimate)
			
			
		xtreg dryrent WNS neighbor_outbreak_count q?precip q?mintemp pop highqual medqual i.year#i.st if year >=1980 & reg ==3, fe cluster(fips)
			return scalar beta3 = _b[WNS]
			return scalar delta3 = _b[neighbor_outbreak_count]
		capture drop newvar3
		gen newvar3 = WNS	+delta3/beta3 * neighbor_outbreak_count

		xtreg agland newvar3   q?precip q?mintemp pop i.year#i.st if year >=1980 & reg ==3, fe cluster(fips)
			return scalar gamma3 = _b[newvar3]
			return scalar alpha3 = _b[_cons]
			lincom _b[newvar3]*delta3/beta3
			return scalar theta3 = r(estimate)			
		
	

	
		restore
	end
	
	simulate beta1 = r(beta1) delta1= r(delta1)   gamma1 = r(gamma1) theta1 = r(theta1) alpha1 = r(alpha1) ///
	beta2 = r(beta2) delta2= r(delta2)   gamma2 = r(gamma2) theta2 = r(theta2) alpha2 = r(alpha2) ///
	beta3 = r(beta3) delta3= r(delta3)   gamma3 = r(gamma3) theta3 = r(theta3) alpha3 = r(alpha3), reps(200) seed(123456789): myboot
 bstat, stat(test) n(257) title("Bootstrap Results--Samples Drawn by County")

restore
save data\tempdata, replace
	
*************************************************************
***************table 8 other lags***********************
***Table of spatial spillovers	
	*****need to run this for each of 10, 25, 50, and 140 km
	****
save data\tempdata, replace
	****counties within XX km
use data\tempdata, clear
xtreg dryrent WNS neighbor_outbreak_count q?precip q?mintemp pop highqual medqual   i.year#i.st if year >=1980 & rg !=6 & rg !=4  & insample ==1, fe cluster(fips)
	outreg2 using output\table8,   keep (WNS neighbor_outbreak_count ) tex replace
	
use data\tempdata, clear
keep if insample ==1
drop if dryrent ==.
keep fips* lndryrent dryrent WNS  q?precip q?mintemp pop highqual medqual year st rg insample 

capture drop _merge
**the  Xkmclean files list the fips codes of each county within X km of a given county
merge m:1 fips using data\140kmclean
*140-up to 99 counties
set more off
forval ii=1/99 {
preserve
keep fips year WNS
rename WNS WNS`ii'
rename fips fipsn`ii'
drop if year ==.
drop if WNS`ii'==.
save data\distance\outbreak_`ii', replace
restore
} 
set more off
forval ii = 1/99 {
capture drop _merge
merge m:1 fipsn`ii' year using data\distance\outbreak_`ii'
drop if fips ==.
replace WNS`ii' = 0 if fipsn`ii'==. & WNS`ii' ==.
}
drop fipsn*
drop _merge

egen WNS_140_n = rowtotal(WNS1	WNS2	WNS3	WNS4	WNS5	WNS6	WNS7	WNS8	WNS9	WNS10 ///
	WNS11	WNS12	WNS13	WNS14	WNS15	WNS16	WNS17	WNS18 ///	
	WNS19	WNS20	WNS21	WNS22	WNS23	WNS24	WNS25	WNS26 ///
	WNS27	WNS28	WNS29	WNS30	WNS31	WNS32	WNS33	WNS34 ///
	WNS35	WNS36	WNS37	WNS38	WNS39	WNS40	WNS41	WNS42 ///
	WNS43	WNS44	WNS45	WNS46	WNS47	WNS48	WNS49	WNS50 ///
	WNS51	WNS52	WNS53	WNS54	WNS55	WNS56	WNS57	WNS58 ///
	WNS59	WNS60	WNS61	WNS62	WNS63	WNS64	WNS65	WNS66 ///
	WNS67	WNS68	WNS69	WNS70	WNS71	WNS72	WNS73	WNS74	///
	WNS75	WNS76	WNS77	WNS78	WNS79	WNS80	WNS81	WNS82	///
	WNS83	WNS84	WNS85	WNS86	WNS87	WNS88	WNS89	WNS90	///
	WNS91	WNS92	WNS93	WNS94	WNS95	WNS96	WNS97	WNS98	///
	WNS99 ///
)
xtreg dryrent WNS WNS_140_n  q?precip q?mintemp pop highqual medqual   i.year#i.st if year >=1980 & rg !=6 & rg !=4 , fe cluster(fips)
    outreg2 using output\table8, keep (WNS WNS_140_n ) tex
 use data\tempdata, clear
 **now for 50km
 **50km has 25 counties
keep if insample ==1
drop if dryrent ==.
keep fips* lndryrent dryrent WNS  q?precip q?mintemp pop highqual medqual year st rg insample 

capture drop _merge
merge m:1 fips using data\50kmclean
set more off
forval ii=1/25 {
preserve
keep fips year WNS
rename WNS WNS`ii'
rename fips fipsn`ii'
drop if year ==.
drop if WNS`ii'==.
save data\distance\outbreak_`ii', replace
restore
} 

set more off
forval ii = 1/25 {
capture drop _merge
merge m:1 fipsn`ii' year using data\distance\outbreak_`ii'
drop if fips ==.
replace WNS`ii' = 0 if fipsn`ii'==. & WNS`ii' ==.
}
drop fipsn*
drop _merge

egen WNS_50_n = rowtotal(	WNS1	WNS2	WNS3	WNS4	WNS5	WNS6	WNS7	WNS8	WNS9	WNS10 ///
	WNS11	WNS12	WNS13	WNS14	WNS15	WNS16	WNS17	WNS18 ///	
	WNS19	WNS20	WNS21	WNS22	WNS23	WNS24	WNS25	 ///
)
set more off
	xtreg dryrent WNS WNS_50_n  q?precip q?mintemp pop highqual medqual   i.year#i.st if year >=1980 & rg !=6 & rg !=4 , fe cluster(fips)
	outreg2 using output\table8,  keep (WNS WNS_50_n ) tex	
	
**25 km has 12 counties	
use data\tempdata, clear	
keep if insample ==1
drop if dryrent ==.
keep fips* lndryrent dryrent WNS  q?precip q?mintemp pop highqual medqual year st rg insample 

capture drop _merge
merge m:1 fips using data\25kmclean

set more off
forval ii=1/12 {
preserve
keep fips year WNS
rename WNS WNS`ii'
rename fips fipsn`ii'
drop if year ==.
drop if WNS`ii'==.
save data\distance\outbreak_`ii', replace
restore
} 
set more off
forval ii = 1/12 {
capture drop _merge
merge m:1 fipsn`ii' year using data\distance\outbreak_`ii'
drop if fips ==.
replace WNS`ii' = 0 if fipsn`ii'==. & WNS`ii' ==.
}
drop fipsn*
drop _merge

egen WNS_25_n = rowtotal(	WNS1	WNS2	WNS3	WNS4	WNS5	WNS6	WNS7	WNS8	WNS9	WNS10 ///
	WNS11	WNS12		 ///
)

xtreg dryrent WNS WNS_25_n  q?precip q?mintemp pop highqual medqual   i.year#i.st if year >=1980 & rg !=6 & rg !=4 , fe cluster(fips)
	outreg2 using output\table8,  keep (WNS WNS_25_n ) tex	

**10 km has 7 counties
	use data\tempdata, clear
	
	keep if insample ==1
drop if dryrent ==.
keep fips* lndryrent dryrent WNS  q?precip q?mintemp pop highqual medqual year st rg insample 

capture drop _merge
merge m:1 fips using data\10kmclean

set more off
forval ii=1/7 {
preserve
keep fips year WNS
rename WNS WNS`ii'
rename fips fipsn`ii'
drop if year ==.
drop if WNS`ii'==.
save data\distance\outbreak_`ii', replace
restore
} 
set more off
forval ii = 1/7 {
capture drop _merge
merge m:1 fipsn`ii' year using data\distance\outbreak_`ii'
drop if fips ==.
replace WNS`ii' = 0 if fipsn`ii'==. & WNS`ii' ==.
}
drop fipsn*
drop _merge

egen WNS_10_n = rowtotal(	WNS1	WNS2	WNS3	WNS4	WNS5	WNS6	WNS7	 ///
			 ///
)

xtreg dryrent WNS WNS_10_n  q?precip q?mintemp pop highqual medqual   i.year#i.st if year >=1980 & rg !=6 & rg !=4 , fe cluster(fips)
	outreg2 using output\table8,  keep (WNS WNS_10_n ) tex
	
******************************************************************8
******************************************************************
use data\tempdata, clear
xtreg dryrent WNS neighbor_outbreak_ind   q?precip q?mintemp pop highqual medqual   i.year#i.st if year >=1980 & rg !=6 & rg !=4 , fe cluster(fips)

gen pop_sample  = pop if dryrent_sample !=.
gen q1precip_sample  = q1precip if dryrent_sample !=.
gen q2precip_sample  = q2precip if dryrent_sample !=.
gen q3precip_sample  = q3precip if dryrent_sample !=.
gen q4precip_sample  = q4precip if dryrent_sample !=.

gen q1mintemp_sample = q1mintemp if dryrent_sample !=.
gen q2mintemp_sample = q2mintemp if dryrent_sample !=.
gen q3mintemp_sample = q3mintemp if dryrent_sample !=.
gen q4mintemp_sample = q4mintemp if dryrent_sample !=.

gen highqual_sample = highqual if dryrent_sample !=.
gen lowqual_sample = lowqual if dryrent_sample !=.
gen medqual_sample = medqual if dryrent_sample !=.
***comparisons for table c1
ttab q1precip_sample q2precip_sample q3precip_sample q4precip_sample ///
q1mintemp_sample q2mintemp_sample q3mintemp_sample q4mintemp_sample ///
pop_sample highqual_sample medqual_sample lowqual_sample ///
  , by(wn)	estout(c(b(fmt(2) star)  se(fmt(3) par)) mlab(,none) collab(,none))	
  
save data\forsumstats, replace  
  bys fips: egen q1mintempavg0 = mean(q1mintemp) if year <=2007 & year >=2005
bys fips: egen q1mintempavg = max(q1mintempavg0)

bys fips: egen q2mintempavg0 = mean(q2mintemp) if year <=2007 & year >=2005
bys fips: egen q2mintempavg = max(q2mintempavg0)

bys fips: egen q3mintempavg0 = mean(q3mintemp) if year <=2007 & year >=2005
bys fips: egen q3mintempavg = max(q3mintempavg0)

bys fips: egen q4mintempavg0 = mean(q4mintemp) if year <=2007 & year >=2005
bys fips: egen q4mintempavg = max(q4mintempavg0)

bys fips: egen q1precipavg0 = mean(q1precip) if year <=2007 & year >=2005
bys fips: egen q1precipavg = max(q1precipavg0)

bys fips: egen q2precipavg0 = mean(q2precip) if year <=2007 & year >=2005
bys fips: egen q2precipavg= max(q2precipavg0)

bys fips: egen q3precipavg0 = mean(q3precip) if year <=2007 & year >=2005
bys fips: egen q3precipavg = max(q3precipavg0)

bys fips: egen q4precipavg0 = mean(q4precip) if year <=2007 & year >=2005
bys fips: egen q4precipavg = max(q4precipavg0)

bys fips: egen popavg0 = mean(pop) if year <=2007 & year >=2005
bys fips: egen popavg= max(popavg0)

bys fips: gen highqual20080 = highqual if year ==2008
bys fips: egen highqual2008 = max(highqual20080)

bys fips: gen medqual20080 = medqual if year ==2008
bys fips: egen medqual2008 = max(medqual20080)

drop q?mintempavg0 q?precipavg0 popavg0 highqual20080 medqual20080


xtreg dryrent WNS  neighbor_outbreak_count    q?precip q?mintemp pop highqual medqual   i.year#i.st if year >=1980 & rg !=6 & rg !=4 & sample==1 , fe cluster(fips)
outreg2 using output\appendix1, keep(WNS neighbor_outbreak_count   ) tex replace


preserve
	collapse (first) wn insample q?mintempavg q?precipavg popavg highqual2008 medqual2008, by(fips)
	ttab q?mintempavg q?precipavg popavg highqual2008 medqual2008 if insample ==1,  by(wn)	estout(c(b(fmt(2) star)  se(fmt(3) par)) mlab(,none) collab(,none))
	tab wn  if insample ==1
	capture drop _weight _pscore
	psmatch2 wn q?mintempavg q?precipavg popavg highqual2008 medqual2008,  noreplacement
	sort fips
	keep fips _weight 
	save data/weights, replace
restore
capture drop _merge  _weight 
merge m:1 fips using data\weights

xtreg dryrent WNS  neighbor_outbreak_count    q?precip q?mintemp pop highqual medqual i.year#i.st if year >=1980 & rg !=6 & rg !=4 [pw = _weight] , fe cluster(fips)
	outreg2 using output\appendix1, keep(WNS neighbor_outbreak_count   ) tex 
	
preserve
	set more off
	collapse (first) wn insample q?mintempavg q?precipavg popavg highqual2008 medqual2008, by(fips)
	ttab q?mintempavg q?precipavg popavg highqual2008 medqual2008 if insample ==1,  by(wn)	estout(c(b(fmt(2) star)  se(fmt(3) par)) mlab(,none) collab(,none))
	tab wn  if insample ==1
	capture drop _weight _pscore
	psmatch2 wn, mahalanobis(q?mintempavg q?precipavg popavg highqual2008 medqual2008) n(5) caliper(2)
	sort fips
	keep fips _weight 
	save data\weights, replace
restore
capture drop _merge  _weight 
merge m:1 fips using data\weights

xtreg dryrent WNS  neighbor_outbreak_count    q?precip q?mintemp pop highqual medqual i.year#i.st if year >=1980 & rg !=6 & rg !=4 [pw = _weight] , fe cluster(fips)
	outreg2 using output\appendix1, keep(WNS neighbor_outbreak_count   ) tex 

	
	preserve
set more off
collapse (first) wn insample q?mintempavg q?precipavg popavg highqual2008 medqual2008, by(fips)
ttab q?mintempavg q?precipavg popavg highqual2008 medqual2008 if insample ==1,  by(wn)	estout(c(b(fmt(2) star)  se(fmt(3) par)) mlab(,none) collab(,none))
tab wn  if insample ==1
capture drop _weight _pscore
psmatch2 wn, mahalanobis(q?mintempavg q?precipavg popavg highqual2008 medqual2008) n(7) caliper(2)
sort fips
keep fips _weight 
save data\weights, replace
restore
capture drop _merge  _weight 

merge m:1 fips using data\weights

xtreg dryrent WNS  neighbor_outbreak_count    q?precip q?mintemp pop highqual medqual i.year#i.st if year >=1980 & rg !=6 & rg !=4 [pw = _weight] , fe cluster(fips)
	outreg2 using output\appendix1, keep(WNS neighbor_outbreak_count   ) tex 


	
preserve
set more off
collapse (first) wn insample q?mintempavg q?precipavg popavg highqual2008 medqual2008, by(fips)
ttab q?mintempavg q?precipavg popavg highqual2008 medqual2008 if insample ==1,  by(wn)	estout(c(b(fmt(2) star)  se(fmt(3) par)) mlab(,none) collab(,none))
tab wn  if insample ==1
capture drop _weight _pscore
psmatch2 wn q?mintempavg q?precipavg popavg highqual2008 medqual2008

sort fips
keep fips _weight _pscore
save data\weights, replace
restore
capture drop _merge  _weight 

merge m:1 fips using data\weights

sum _pscore, det
	
	set more off
	xtreg dryrent WNS neighbor_outbreak_count     q?precip q?mintemp pop highqual medqual i.year#i.st if year >=1980 & rg !=6 & rg !=4 &_pscore>=.005& _pscore<=.22 , fe cluster(fips)
	outreg2 using output\appendix1, keep(WNS neighbor_outbreak_count   ) tex 	

***************************************************************************

******************************************************************************
******summary stats, table 1****
use data/forsumstats, clear
eststo sumvars: estpost summarize dryrent_sample agland_sample  corny_sample wheaty_sample soyy_sample exp_land_sample  q?precip_sample q?mintemp_sample pop_sample highqual_sample medqual_sample lowqual_sample WNS_sample WNS_n_count_sample WNS_n_n_sample
	esttab sumvars  using output\sumstats.tex, replace cell("mean(fmt(2)) sd(fmt(2)) count(fmt(0))") 

***************************************************************************
***with lagged weather for the appendix table c3
	sort fips year
	xtreg dryrent WNS  neighbor_outbreak_count    q?precip q?mintemp l.q?precip l.q?mintemp pop highqual medqual   i.year#i.st if year >=1980 & rg !=6 & rg !=4 , fe cluster(fips)
	outreg2 using output\appendix2, replace keep (WNS  neighbor_outbreak_count   )  tex
	xtreg dryrent WNS neighbor_outbreak_ind     q?precip q?mintemp l.q?precip l.q?mintemp  pop highqual medqual  i.year#i.st if year >=1980 & rg !=6 & rg !=4 , fe cluster(fips)
	outreg2 using output\appendix2,  keep (WNS  neighbor_outbreak_count neighbor_outbreak_ind   )   tex
	xtreg lndryrent WNS neighbor_outbreak_count     q?precip q?mintemp l.q?precip l.q?mintemp  pop highqual medqual  i.year#i.st if year >=1980 & rg !=6 & rg !=4 , fe cluster(fips)
	outreg2 using output\appendix2,  keep (WNS neighbor_outbreak_count   )   tex
xtreg dryrent WNS      q?precip q?mintemp pop l.q?precip l.q?mintemp  highqual medqual  i.year#i.st if year >=1980 & rg !=6 & rg !=4 , fe cluster(fips)
	outreg2 using output\appendix2,  keep (WNS  neighbor_outbreak_count   )   tex
******************************************************************************
****All LCC, appendix table c4

xtreg dryrent WNS  neighbor_outbreak_count    q?precip q?mintemp pop qual*   i.year#i.st if year >=1980 & rg !=6 & rg !=4 & sample==1 , fe cluster(fips)
	outreg2 using output\appendix3, replace keep (WNS neighbor_outbreak_count    ) tex
	xtreg dryrent WNS  neighbor_outbreak_ind    q?precip q?mintemp pop  qual*     i.year#i.st if year >=1980 & rg !=6 & rg !=4  & sample==1, fe cluster(fips)
	outreg2 using output\appendix3,  keep (WNS neighbor_outbreak_count  neighbor_outbreak_ind   ) tex
	xtreg lndryrent WNS  neighbor_outbreak_count    q?precip q?mintemp pop  qual*     i.year#i.st if year >=1980 & rg !=6 & rg !=4  & sample==1, fe cluster(fips)
	outreg2 using output\appendix3,  keep (WNS neighbor_outbreak_count   ) tex
xtreg dryrent WNS   q?precip q?mintemp pop qual*      i.year#i.st if year >=1980 & rg !=6 & rg !=4  & sample==1, fe cluster(fips)
	outreg2 using output\appendix3,  keep (WNS neighbor_outbreak_count   ) tex

************************************************************
**************************************************************
*********************************************************
*****************SIMULATIONS****************************************
****output used in appendix tables D1 and D2, as well as Figures 4 and 6

set more off
	***price change
	xtreg dryrent WNS neighbor_outbreak_count q?precip q?mintemp pop i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1 ,  fe  cluster(fips)
			scalar beta = _b[WNS]
			scalar delta = _b[neighbor_outbreak_count]
			capture drop newvar
		gen newvar = WNS	+delta/beta * neighbor_outbreak_count
	xtreg agland  newvar q?precip q?mintemp pop i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1,  fe cluster(fips) 
		scalar gamma = _b[newvar]
		scalar alpha = _b[_cons]
		lincom _b[newvar]*delta/beta
		scalar theta = r(estimate)	

	**base amounts
	sum agland if year ==2007 & insample ==1
	scalar land_mean = r(mean)
	display land_mean
	sum dryrent if year==2008 & insample ==1
	scalar baserent = r(mean)
	sum agland if year ==2007 
	scalar land_all = r(mean)
	
	**parameters for direct impact
	scalar p0 = baserent 
	scalar p1_direct = baserent + beta
	scalar p1_spill1 = baserent + 1*delta
	scalar p1_spill2 = baserent + 2*delta
	scalar p1_spill3 = baserent + 3*delta
	scalar p1_spill4 = baserent + 4*delta
	scalar p1_spill5 = baserent + 5*delta
	scalar q0 = land_mean
	scalar q1_direct = land_mean + gamma
	scalar q1_spill1 = land_mean + 1*theta
	scalar q1_spill2 = land_mean + 2*theta
	scalar q1_spill3 = land_mean + 3*theta
	scalar q1_spill4 = land_mean + 4*theta
	scalar q1_spill5 = land_mean + 5*theta
	***********CHANGE eps HERE:
	scalar eps = .45
	scalar b = p0/q0*1/eps
	disp b
	
	
	scalar d_direct = (p1_direct-p0)/(q1_direct-q0)
	scalar d_spill1 = (p1_spill1-p0)/(q1_spill1-q0)
	scalar d_spill2 = (p1_spill2-p0)/(q1_spill2-q0)
	disp d_direct
	disp d_spill1
	disp d_spill2
	scalar eta_direct = (1/d_direct) * p0/q0
	scalar eta_spill1 = (1/d_spill1) * p0/q0
	scalar eta_spill2 = (1/d_spill2) * p0/q0
	disp eta_spill1
	disp eta_spill2
	scalar c_direct = p0 - d_direct*q0
	scalar c_spill1 = p0 - d_spill1*q0
	scalar c_spill2 = p0 - d_spill2*q0
	scalar c_spill3 = p0 - d_spill2*q0
	scalar c_spill4 = p0 - d_spill2*q0
	scalar c_spill5 = p0 - d_spill2*q0
	disp c_direct
	disp c_spill1
	disp c_spill2
	disp c_spill3
	disp c_spill4
	disp c_spill5
	
	scalar a = p0+b*q0
	disp a
	scalar a2 = p0*(1+1/eps)
	disp a2	
	disp p0 
	disp p1_direct
	disp p1_spill1
	disp q0
	disp q1_direct
	disp q1_spill1
	disp q1_spill2
	***impact
	
	scalar K_par_direct = p0 - p1_direct + p0/q0 * 1/eps * (q0 - q1_direct)
	scalar K_par_spill1 = p0 - p1_spill1 + p0/q0 * 1/eps * (q0 - q1_spill1)
	scalar K_par_spill2 = p0 - p1_spill2 + p0/q0 * 1/eps * (q0 - q1_spill2)
	scalar K_par_spill3 = p0 - p1_spill3 + p0/q0 * 1/eps * (q0 - q1_spill3)
	scalar K_par_spill4 = p0 - p1_spill4 + p0/q0 * 1/eps * (q0 - q1_spill4)
	scalar K_par_spill5 = p0 - p1_spill5 + p0/q0 * 1/eps * (q0 - q1_spill5)
	
	**SET CHANGE IN SLOPE HERE
	**parallel:
	scalar db_direct = 0
	scalar db_spill1 = 0
	scalar db_spill2 = 0
	scalar db_spill3 = 0
	scalar db_spill4 = 0
	scalar db_spill5 = 0
		
	**to do a pivot around y-intercept:
	*scalar db = K_par/q1
	/*
	scalar db_direct = K_par_direct/q1_direct	
	scalar db_spill1 = K_par_spill1/q1_spill1	
	scalar db_spill2 = K_par_spill2/q1_spill2	
	scalar db_spill3 = K_par_spill3/q1_spill3	
	scalar db_spill4 = K_par_spill4/q1_spill4	
	scalar db_spill5 = K_par_spill5/q1_spill5	
	*/
	scalar K_direct = p0 - p1_direct + p0/q0 * 1/eps * (q0 - q1_direct) - db_direct * q1_direct
	
	scalar K_spill1 = p0 - p1_spill1 + p0/q0 * 1/eps * (q0 - q1_spill1) - db_spill1 * q1_spill1
	scalar K_spill2 = p0 - p1_spill2 + p0/q0 * 1/eps * (q0 - q1_spill2) - db_spill2 * q1_spill2
	scalar K_spill3 = p0 - p1_spill3 + p0/q0 * 1/eps * (q0 - q1_spill3) - db_spill3 * q1_spill3
	scalar K_spill4 = p0 - p1_spill4 + p0/q0 * 1/eps * (q0 - q1_spill4) - db_spill4 * q1_spill4
	scalar K_spill5 = p0 - p1_spill5 + p0/q0 * 1/eps * (q0 - q1_spill5) - db_spill5 * q1_spill5
	
	scalar Kcheck_direct = p0 - p1_direct + p0/eps- p0/eps *q1_direct/q0	
	scalar Kcheck_spill1 = p0 - p1_spill1 + p0/eps- p0/eps *q1_spill1/q0	
	scalar Kcheck_spill2 = p0 - p1_spill2 + p0/eps- p0/eps *q1_spill2/q0	
	
	**to double (comment out for other sims):
	/*
	scalar K_direct = 2*K_par_direct
	scalar K_spill1 = 2*K_par_spill1
	scalar K_spill2 = 2*K_par_spill2
	scalar K_spill3 = 2*K_par_spill3
	scalar K_spill4 = 2*K_par_spill4
	scalar K_spill5 = 2*K_par_spill5
	*/
	**pivot around new equilibrium:
	*scalar K = a-p1
*	scalar K = 0 
	
	disp  K_direct
	disp K_spill1
	disp K_spill2
	disp K_spill3
	disp K_spill4
	disp K_spill5
	**K appears to scale linearly
	disp  Kcheck_direct
	disp  Kcheck_spill1
	disp  Kcheck_spill2
	disp db_direct
	
	**implied db from K
	*intercept:
	scalar int_test_direct = a-K_direct
	scalar slope_test_direct = (a - K_direct - p1_direct)/q1_direct
	scalar db_test_direct = b - slope_test_direct
	
	scalar int_test_spill1 = a-K_spill1
	scalar slope_test_spill1 = (a - K_spill1 - p1_spill1)/q1_spill1
	scalar db_test_spill1 = b - slope_test_spill1
	
	scalar int_test_spill2 = a-K_spill2
	scalar slope_test_spill2 = (a - K_spill2 - p1_spill2)/q1_spill2
	scalar db_test_spill2 = b - slope_test_spill2
	
	disp int_test_direct
	disp slope_test_direct
	disp db_test_direct
	
	disp int_test_spill1
	disp slope_test_spill1
	disp db_test_spill1
	
	disp int_test_spill2
	disp slope_test_spill2
	disp db_test_spill2
	
		
	**need to account for c in this calculation
	scalar Loss_direct = 0.5*(a-c_direct)*q0-0.5*(a-K_direct-c_direct)*q1_direct
	scalar Loss_p_direct = q1_direct*K_direct+0.5*(q0-q1_direct)*K_direct
	disp Loss_direct 
	disp Loss_p_direct
	
	scalar Loss_spill1 = 0.5*(a-c_spill1)*q0-0.5*(a-K_spill1-c_spill1)*q1_spill1
	scalar Loss_spill2 = 0.5*(a-c_spill2)*q0-0.5*(a-K_spill2-c_spill2)*q1_spill2
	scalar Loss_spill3 = 0.5*(a-c_spill3)*q0-0.5*(a-K_spill3-c_spill3)*q1_spill3
	scalar Loss_spill4 = 0.5*(a-c_spill4)*q0-0.5*(a-K_spill4-c_spill4)*q1_spill4
	scalar Loss_spill5 = 0.5*(a-c_spill5)*q0-0.5*(a-K_spill5-c_spill5)*q1_spill5
	di %20.0g scalar(Loss_spill1)
	di %20.0g scalar(Loss_spill2)

	disp Loss_spill1
	disp Loss_spill2
	disp Loss_spill3
	disp Loss_spill4
	disp Loss_spill5
			
	scalar qbar_direct = -c_direct/d_direct
	disp qbar_direct
	scalar perzero_direct = qbar_direct/q0
	disp perzero_direct
	
	scalar qbar_spill = -c_spill1/d_spill1
	disp qbar_spill
	scalar perzero_spill = qbar_spill/q0
	disp perzero_spill

	
	*percent of base welfare:
	scalar WF_0_direct = .5*(a-c_direct)*q0-.5*(-c_direct)*qbar_direct
	disp WF_0_direct
	scalar WF_0 = .5*(a-c_direct)*q0-.5*(-c_direct)*qbar_direct
	scalar WF_per_base_direct = Loss_direct/WF_0_direct
	disp WF_per_base_direct
	
	scalar WF_0_spill1 = .5*(a-c_spill1)*q0-.5*(-c_spill1)*qbar_spill
	disp WF_0_spill1
	scalar WF_per_base_spill1 = Loss_spill1/WF_0_spill1
	disp WF_per_base_spill1
	
	*per base acre
	scalar lossperacre_direct = Loss_direct/q0
	disp lossperacre_direct 

	scalar lossperacre_spill1 = Loss_spill1/q0
	disp lossperacre_spill1 
	
	scalar loss_direct_check = Loss_direct*178
	disp loss_direct_check
	
***Total direct impact
*as of end of data, max(WNS), gen direct cost = Loss_direct * that var, then summed
preserve
keep if insample==1
bys fips: egen direct_WNS = max(WNS)
collapse (first) direct_WNS insample, by (fips )
sum direct_WNS, det
gen direct_cost_co = direct_WNS*Loss_direct
egen direct_cost_total = sum(direct_cost_co)
format %20.10f direct_cost_total
sum direct_cost_total direct_WNS
restore

****Total indirect impact
**as of end of data	

preserve
keep if insample==1
bys fips: egen spill_WNS = max(neighbor_outbreak_count)
collapse (first) spill_WNS insample, by (fips )
sum spill_WNS, det
tab spill_WNS
gen spill_cost_co = Loss_spill1 if spill_WNS == 1
replace spill_cost_co = Loss_spill2 if spill_WNS == 2
replace spill_cost_co = Loss_spill3 if spill_WNS == 3
replace spill_cost_co = Loss_spill4 if spill_WNS == 4
replace spill_cost_co = Loss_spill5 if spill_WNS == 5
egen spill_cost_total = sum(spill_cost_co)
format %20.10f spill_cost_total 
sum spill_cost_total spill_WNS
restore

preserve
	matrix test = (K_direct, K_spill1 ,K_spill2 ,K_spill3 ,K_spill4 ,K_spill5 , Loss_direct, Loss_spill1, Loss_spill2,Loss_spill3,Loss_spill4,Loss_spill5, WF_per_base_direct, WF_per_base_spill1, WF_0 )
matrix list test
*****change eps in 2 places	
	**with bootstrap***
	capture program drop myboot
		 program define myboot, rclass
		 preserve
		 bsample, cluster(fips)
		 set more off
		***price change
	xtreg dryrent WNS neighbor_outbreak_count q?precip q?mintemp pop i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1 ,  fe  cluster(fips)
			scalar beta = _b[WNS]
			scalar delta = _b[neighbor_outbreak_count]
			capture drop newvar
		gen newvar = WNS	+delta/beta * neighbor_outbreak_count
	xtreg agland  newvar q?precip q?mintemp pop i.year#i.st if year >=1980 & rg !=6 & rg !=4 & insample ==1,  fe cluster(fips) 
		scalar gamma = _b[newvar]
		scalar alpha = _b[_cons]
		lincom _b[newvar]*delta/beta
		scalar theta = r(estimate)	

		sum agland if year ==2007 & insample ==1
		scalar land_mean = r(mean)
		sum dryrent if year==2008 & insample ==1
		scalar baserent = r(mean)
		sum agland if year ==2007 
		scalar land_all = r(mean)
		
		**parameters
		scalar p0 = baserent 
		scalar q0 = land_mean
		scalar p1_direct = baserent + beta
		scalar p1_spill1 = baserent + 1*delta
		scalar p1_spill2 = baserent + 2*delta
		scalar p1_spill3 = baserent + 3*delta
		scalar p1_spill4 = baserent + 4*delta
		scalar p1_spill5 = baserent + 5*delta
		
		scalar q1_direct = land_mean + gamma
		scalar q1_spill1 = land_mean + 1*theta
		scalar q1_spill2 = land_mean + 2*theta
		scalar q1_spill3 = land_mean + 3*theta
		scalar q1_spill4 = land_mean + 4*theta
		scalar q1_spill5 = land_mean + 5*theta
			
		/*
		scalar db = 0
		scalar eta = (p0/q0)*((q1-q0)/(p1-p0))
		scalar c =  p0*(1-1/eta) 
	
		scalar a = p0*(1+1/eps)
		scalar d = (1/eta)*p0/q0
		scalar qbar = -c/d
		scalar eps = .55
		scalar b = p0/q0*1/eps
		scalar WF_base = .5*((1+1/eps)-(1-1/eta))*p0*q0-.5*(-c)*qbar
		scalar K0 = p0 - p1
		
		return scalar K = p0 - p1 + p0/q0 * 1/eps * (q0 - q1) + db * q1
		return scalar Loss = q1*K+0.5*(q0-q1)*K
		return scalar Loss0 = q1*K0+0.5*(q0-q1)*K0
		return scalar WF_per_base = Loss/WF_base
		*/
		scalar d = (p1_direct-p0)/(q1_direct-q0)
		scalar c = p0 - d*q0 
		scalar eps = .45
		scalar b = p0/q0*1/eps
		scalar a = p0+b*q0
		
		scalar qbar = -c/d
		/*
		scalar K_par_direct = p0 - p1_direct + p0/q0 * 1/eps * (q0 - q1_direct)
		scalar K_par_spill1 = p0 - p1_spill1 + p0/q0 * 1/eps * (q0 - q1_spill1)
		scalar K_par_spill2 = p0 - p1_spill2 + p0/q0 * 1/eps * (q0 - q1_spill2)
		scalar K_par_spill3 = p0 - p1_spill3 + p0/q0 * 1/eps * (q0 - q1_spill3)
		scalar K_par_spill4 = p0 - p1_spill4 + p0/q0 * 1/eps * (q0 - q1_spill4)
		scalar K_par_spill5 = p0 - p1_spill5 + p0/q0 * 1/eps * (q0 - q1_spill5)
		*/
		return scalar WF_0 = .5*(a-c)*q0-.5*(-c)*qbar
		
		/*
		scalar db_direct = K_par_direct/q1_direct
		scalar db_spill1 = K_par_spill1/q1_spill1
		scalar db_spill2 = K_par_spill2/q1_spill2	
		scalar db_spill3 = K_par_spill3/q1_spill3	
		scalar db_spill4 = K_par_spill4/q1_spill4	
		scalar db_spill5 = K_par_spill5/q1_spill5
		*/
		
		scalar db_direct = 0
		scalar db_spill1 = 0
		scalar db_spill2 = 0
		scalar db_spill3 = 0
		scalar db_spill4 = 0
		scalar db_spill5 = 0
		
		return scalar K_direct = p0 - p1_direct + p0/q0 * 1/eps * (q0 - q1_direct) - db_direct * q1_direct
		return scalar K_spill1 = p0 - p1_spill1 + p0/q0 * 1/eps * (q0 - q1_spill1) - db_spill1 * q1_spill1
		return scalar K_spill2 = p0 - p1_spill2 + p0/q0 * 1/eps * (q0 - q1_spill2) - db_spill2 * q1_spill2
		return scalar K_spill3 = p0 - p1_spill3 + p0/q0 * 1/eps * (q0 - q1_spill3) - db_spill3 * q1_spill3
		return scalar K_spill4 = p0 - p1_spill4 + p0/q0 * 1/eps * (q0 - q1_spill4) - db_spill4 * q1_spill4
		return scalar K_spill5 = p0 - p1_spill5 + p0/q0 * 1/eps * (q0 - q1_spill5) - db_spill5 * q1_spill5

		/*
		return scalar K_direct = 2*K_par_direct
		return scalar K_spill1 = 2*K_par_spill1
		return scalar K_spill2 = 2*K_par_spill2
		return scalar K_spill3 = 2*K_par_spill3
		return scalar K_spill4 = 2*K_par_spill4
		return scalar K_spill5 = 2*K_par_spill5
*/
	    return scalar Loss_direct = 0.5*(a-c)*q0-0.5*(a-K_direct-c)*q1_direct
		return scalar Loss_spill1 = 0.5*(a-c)*q0-0.5*(a-K_spill1-c)*q1_spill1
		return scalar Loss_spill2 = 0.5*(a-c)*q0-0.5*(a-K_spill2-c)*q1_spill2
		return scalar Loss_spill3 = 0.5*(a-c)*q0-0.5*(a-K_spill3-c)*q1_spill3
		return scalar Loss_spill4 = 0.5*(a-c)*q0-0.5*(a-K_spill4-c)*q1_spill4
		return scalar Loss_spill5 = 0.5*(a-c)*q0-0.5*(a-K_spill5-c)*q1_spill5
	    return scalar WF_per_base_direct = Loss_direct/WF_0
		return scalar WF_per_base_spill1 = Loss_spill1/WF_0
		
		restore
	end
	
	simulate K_direct = r(K_direct) K_spill1 = r(K_spill1) K_spill2 = r(K_spill2) ///
	K_spill3 = r(K_spill3) K_spill4 = r(K_spill4) K_spill5 = r(K_spill5) Loss_direct = r(Loss_direct) Loss_spill1 = r(Loss_spill1) ///
	Loss_spill2 = r(Loss_spill2) Loss_spill3 = r(Loss_spill3) Loss_spill4 = r(Loss_spill4) Loss_spill5 = r(Loss_spill5) ///
	WF_per_base_direct = r(WF_per_base_direct) WF_per_base_spill1 = r(WF_per_base_spill1) WF_0 = r(WF_0) , reps(100) seed(123456789): myboot
 bstat, stat(test) n(257) title("Bootstrap Results--Samples Drawn by County")
restore
	

***********************************
***figure 3
use data/bat_dataset, clear
bys fips: egen maxtimeafter = max(event2)

bys fips: egen obs = count(event)
capture drop obs2
bys fips: egen obs2 = count(dryrent) if obs !=. & event >=-3 & event <=3 & event !=.	

	reg dryrent b42.event2 if year >=1980 & rg !=4   & rg !=6    & wn ==1 & event >= -5 & event<=5 ,  cluster(fips) 
	est store eve2 
	coefplot eve2, keep(*.event2 ) vertical ci(95) coeflabels(30.event2= "-12" 31.event2 ="-11" 32.event2 ="-10" 33.event2 ="-9" ///
	34.event2 ="-8" 35.event2 ="-7" 36.event2 ="-6" 37.event2 ="-5" 38.event2 ="-4" 39.event2 ="-3" 40.event2 ="-2" ///
	41.event2 ="-1" 43.event2 ="1" 44.event2 ="2" 45.event2 ="3" 46.event2 ="4" 47.event2 ="5" 48.event2 ="6"  ///
	49.event2 ="7" 50.event2 ="8") xtitle("Years Since WNS Arrival") ytitle("Change in Cash Rent per Acre")
		
***************************************************************************888

******************************************************************************
********for the map (figure 2:
preserve
	xtreg dryrent WNS neighbor_outbreak_ind q?precip q?mintemp pop i.year#i.st if year >=1980 & rg !=6 & rg !=4 , fe cluster(fips)
	gen sample = e(sample)
	*bys fips: egen insample = max(sample)
	keep if sample ==1
	bys fips: gen arrive = year if wn ==1 & sample==1 & WNS==1
	bys fips: egen arrival = min(arrive)
	gen treatment = 0 if arrival ==.
	replace treatment = 1 if arrival !=.
	replace treatment = 0 if arrival ==.
	rename county subregion
	capture drop region
	rename state region
	
	collapse (first)  region subregion arrival treatment, by(fips)
	rename arrival year
	
	bys fips year: gen test = _n
	sum test
	*rename fips geoid
	gen geoid = string(fips,"%05.0f")
	*gen geoid = fips
	keep geoid year region subregion treatment
	tab treatment
	tab year
	save data\tomap.dta, replace
restore		
**then run R code, ggplot_sf_map_bats
* 178 have wns in the sample--176 have an arrival date during the sample


******end of replication code in stata************8

